Symplectic algorithms for simulations of rigid body systems using the exact solution 

of free motion 



o 
o 

(N 



o 

CO 

-4-* 



C 

o 

(N 
> 
^- 
O 
^t- 
(N 

o 

-I— > 



i 

c 

o 
o 



X 
S3 



Ramses van Zon and Jeremy Schofield 
Chemical Physics Theory Group, Department of Chemistry, University of Toronto, 
80 Saint George Street, Toronto, Ontario, Canada M5S 3H6 
(Dated: May 17, 2007) 

Elegant integration schemes ol second and fourth order for simulations of rigid body systems are 
presented which treat translational and rotational motion on the same footing. This is made possible 
by a recent implementation of the exact solution of free rigid body motion. The two schemes are 
time-reversible, symplectic, and exactly respect conservation principles for both the total linear and 
angular momentum vectors. Simulations of simple test systems show that the second order scheme 
is stable and conserves all constants of the motion to high precision. Furthermore, the schemes are 
demonstrated to be more accurate and efficient than existing methods, except for high densities, in 
which case the second order scheme performs at least as well, showing their general applicability. 
Finally, it is demonstrated that the fourth order scheme is more efficient than the second order 
scheme provided the time step is smaller than a system-dependent threshold value. 

PACS numbers: 45.10.-b, 02.70.Ns, 45.40.-f, 61.20.Ja 



I. INTRODUCTION 

Systems of rigid bodies can be used to model phenom- 
ena on a wide ranee of length scales, from the dynam- 
ics of molecules [l], 0, ID, H[, polymers and other complex 
systems|5], to robotics on a macroscopic level[6|. On even 
larger scales, many astrophysical objects such as plan- 
ets, satellites, and space crafts can be regarded as rigid 
bodies @, § • Much recent work has been devoted to devis- 
ing time-reversible and symplectic numerical integrators 
for rigid systems [l|, 0> S, H|. I Q this paper, we present 
second and fourth order symplectic integration schemes 
for general interacting rigid bodies that make use of a 
recent numerical implementation of the exact solution of 
free rig id bodies [3. The schemes do not require speci- 
fying orientational parameters such as Euler angles, nor 
extending the state space, nor the use of quaternions and 
constraint conditions, and have a symplectic structure on 
the standard phase space of rigid systems. Besides being 
aesthetically appealing, the schemes are more efficient at 
a given level of accuracy in many cases, and never less 
efficient. Thus they are good multi-purpose integrators. 

The paper is structured as follows. In Sec. [Til the dy- 
namics of rigid body systems is briefly reviewed. The 
algorithms that apply to such systems are proposed in 
Sec. lIIIl while their symplecticity is demonstrated in 
Sec. lIVI ScctionfVl contains numerical tests on three 
model systems: a free asymmetric rotating body (|V A[) . 
a water molecule in an electric field (|VB[) and a model of 
liquid water (|V Cp . The paper ends with the conclusions 
in Seclvn 



II. DYNAMICS OF RIGID BODY SYSTEMS 

Consider a system consisting of N rigid bodies, where 
each body i has a mass and is described by its center- 
of-mass position and its attitude matrix . The rows 



of Ai are the directions of the principle axes of the body in 
the lab frame. The position of a point a with coordinates 
f a in the principal axis frame of body i is given by r a — 
qi + Af r a . In addition, each body has momentum 



p, ; = niiCii 



(1) 



and angular velocity u>i. In index notation, uii is related 
to the time derivative of the attitude matrix vial 10] 



^ £bac(Ui)a = ^ (hi)ab(Ai) ac 



(2) 



a—x,y,z 



a—x,y,z 



where ebac is the Levi-Civita symbol[n|. A related quan- 
tity is the angular momentum = \iU>i, whose impor- 
tance lies in the fact that the total angular momentum 
Lt = x Pi + -^i) is a conserved quantity in rota- 

tionally symmetric systems. Here, li is the moment of 
inertia tensor lj = A^I^A^, where 1^ = di&g(Ii X , I{ y , I{ Z ) 
and Ii x , Ii y and Ii Z are the principal moments of inertia. 
The dynamics of the system is given bv[ll| 



Pi = Fi, 

Li = Ti, 



(3a) 
(3b) 



where the forces F^ and torques t; are assumed to be 
functions of {q^} and {Aj} only. The forces may derive 
from an interaction potential V between sites, such that 
the force on site a of body i is Fi a = — d ria V. These 
forces then determine the body forces and torques: 



Fi ^ Fi Q , 

Q = l 

Ti = ^2{r ia -qi)x Fi 



(4a) 
(4b) 



where rii is the number of interaction sites on body i. 
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III. INTEGRATION ALGORITHMS 

The proposed algorithm to integrate Eqs. (fT|)- ([3b"l) is 
based on the exact solution of torque-free dynamics and 
will be presented first, after which its derivation will be 
given. 

In the algorithm, each body i is specified by the set 
(qi, Pi, A,-,Lj), and two different evolution operators are 
defined. The first is the exact free evolution operator 
ipl over a time h, under which Lj and Pi are invariant 
while and Ai evolve according to: 



(o) 



q, + hpi/nti, 



(5a) 
(5b) 



where the matrix Pi(/i) (specified below) propagates Ai 
from time t to t* =t + h. The second evolution operator, 
denoted by , evolves the momenta due to interactions: 



(i) 



Pi 

L; 



flTi 



(6a) 
(6b) 



while qi and Ai remain unchanged. Before applying this 
operator, the forces and torques must be computed using 
standard techniques [l2l Il3j . 

In the integration schemes, the evolution over a time 
t is replaced by evolution over t/h time steps of dura- 
tion h. In each time step, the two evolution operators 
(/A ) and are applied in sequence. For the second 
order integrator, this sequence is that of a generalized 
Verlet scheme, i.e. 



(7) 



This second order symplectic integration scheme using 
the exact rotation matrix will be referred to as SIER2. 

For the fourth order variant of the integrator, de- 
noted as SIER4, the extended Forest-Ruth-like scheme 
of Omelyan et a?. [13] is used, in which the true evolution 
is replaced by 

,„(0) ,.(1) .JO) (1) (0) (1) (0) (1) n(h 5 s ( os 

Vht ^h 2 ^h 3 <Ph 4 ^hs Vhi Vh 3 H> hl +U(h ), (8) 



with hi = j — /12, he, = h — 2(h\ + h%), h\ = 
0.1720865590295143/i, h 2 = 0.5915620307551568 h, and 
h 3 = -0.1616217622107222/j. 

The most novel aspect of these integrators is the use 
of the exact solution of Pi in y>(°) [l5|. Its form is [Tol] 



P(h) = R^L^^^R^L), 



(9) 



where the i dependence has been dropped with the under- 
standing that each body i has its own matrix Pj to prop- 
agate its specific attitude matrix Ai . In Eq. © , L = AL 
is the angular momentum in the principal axes frame at 
time t, L* is the same quantity at time t*, and Ri is 



Ri(L) 



1 



L X L Z —LLy L±L X 



LL 



LyL z 

-Li 



-LL X 




L±Ly 

L i L, 



with L± = (Ll+L^) 1 / 2 and L = |L|. While L is known at 

the start of the time step, L* can be computed from the 
free solution of the Euler equations^ [icil]. which involve 
Jacobi elliptic functions. Fortunately, efficient routines 
are available to evaluate such functions [l6j|. 

The matrix R2NA in Eq. ((9|) is a rotation around the 
z-axis by an angle [10] 



ip = ah + ac I w(s) ds, 



(10) 



with a = L/I z , w(s) = L 2 /L\{s), c = 2I Z E/L 2 - 1, 
where E = L • r x L/2. 

The explicit solution of the integral in Eq. (fT0|) re- 
quires the evaluation of elliptic integrals and theta 
functions which can become somewhat of a com- 
putational burden [l7|. Instead, standard numerical ap- 
proaches could be used, such as quadrature-based meth- 
ods. However, quadrature-based methods involve evalu- 
ating additional elliptic functions at several intermediate 
points in the interval (t, t*). The following method of 
numerically computing the integral does not suffer from 
this drawback. 

For small time steps h, the integral in Eq. (fit)]) can 
be expressed using the integral of the polynomial ap- 
proximation of the integrand uu(s) found by Hermite 
interpolation[l8j], yielding 



w(s) ds = y — ^ — „ [w {i x > + (-y 

j!(n-j)!(2n)! L 



(11) 

plus corrections of 0(h 2n+1 ). Here, = d j w(t)/dV, 



and 



w^(t M ). In the present case, 10^ and w 



can be expressed in terms of L and , respectively, using 
the Euler equations recursively. To be more precise, 



2{I X Iy)L x L y L z 



IxIyL 



while in general 



w 



W) = Sj(w) x 



OF 

w^- 1 ' en- 



it j is even, 
if j is odd, 



where the Sj(w) are polynomials in w. The first few are 
So = w, Si = 1, and S2 = 2a + 4(6c — a)w + 6(c — 
b)cw 2 — 8c 2 iu 3 , where a = (I z /I y — 1)(I Z /I X — 1) and 
b = 2 — I z /I y — I z /Ix- The Sj satisfy the recursion Sj = 
Sj^idhj / dw + 2hjdSj-i/dw with hj — 2u>(l — w)(a + 
bcw + c 2 w 2 ) for even j, and Sj = dSj-x/dw for odd j[19|. 

Since odd time derivatives change sign under time re- 
versal, each term in Eq. (jlip is invariant under time rever- 
sal, so this approximation preserves the time-reversibility 
of the integrators. By increasing n until ip has converged, 
one obtains a "numerically exact" result. 

It should be remarked that the equations above apply 
if 2E > I y L 2 . Otherwise, one must perform a rotation 
U* in the principal axes frame to exchange the x and z 
direction and reverse the y direction [lOj]. 
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IV. SYMPLECTICITY 

It is clear that Eqs. (|5b[) — (fT|) generalize the Verlet inte- 
gration algorithm to rigid systems. But it is not obvious 
that this integrator [or the higher-order variant in Eq. (JHJ] 
is symplectic, since symplccticity is a property of the 
equations of motion in Hamiltonian mechanics. In this 
formulationfll], the state of the system is described by 
generalized coordinates and momenta conjugate to these 
coordinates. The structure of the equations of motion 
is symplectic if T = JdrTi., where T is the phase point 
consisting of the generalized coordinates and momenta, 
and Ti is the Hamiltonian ll|. Here J is defined to be 
J = J), where 1 is a 67V x QN unit matrix. For 
time-independent Ti, the evolution over a time h may 
be written as T(h) = exp[Lh]T (0) , where the Liouvillian 
L = —dr"H ■ J<9r- The operator expjX/i] is called the 
propagator and is symplectic due to the form of L. 

While there are several ways to parametrize the at- 
titude matrix A, three generalized coordinates are al- 
ways required, denoted here by ■di = ($u , da, 'da). 
From the Lagrangian of interacting rigid bodies C = 
3 LiIi( TO dqi| 2 +Ui ■ k(Di) - V({qj,0j}), one finds 



H 



It. 



7Ti • M, 1 7T ; 



VdqjM), (12) 



where the generalized momenta are given by pi = m^q 
and TTi = Mji?j. Here, Mj is a i?i-dependent ma- 
trix given by Mj = NiljN^, where, in turn, N a {, = 
\ebcdA ec dA e d/ 'd'&a- To derive the integration schemes 
(0) and ijHJ), the Hamiltonian is split into Ti^ = 
Vdqjrfj}) and = H - Ti {1) . From the H u \ the 
partial Liouvillians = —dpTi^ ■ Jc?r and symplectic 
propagators exp[L^'h] are constructed. The propagator 
exp[i> )/i] evolves the rigid system freely over a time h 
and thus coincides with the propagator ip^ in Eq. {5bJ. 
The action of the propagator exp[L^h] on a phase point 
r can also be determined: 



(qi,<?i,Pi,7r<) = (qi,i?i,Pi - hdq t V,m - hd^V), 



from which it follows that [cf. Eq. ([FBI 



e LWh U = U 



hN^d^V = Li + hTi =^ 1} L 



Hence both ip^ and ip^ are symplectic operators. Since 
a succession of symplectic operators is also symplectic, 
the SIER integrators in Eqs. ([7]) and ([8]) are symplec- 
tic. It is straightforward to show that SIER2 and SIER4 
schemes are second and fourth order, respectively (l2l . 

Both in Eq. (|5bj) and <£>j, in Eq. (|6bj) strictly con- 
serve the angular momentum vector for spherically 
symmetric systems. The conservation of Tit under ip^ 
arises because Lj is a constant of the exact free motion, 
whereas the equality ipy'Lix = Lt is a consequence of the 



rotational invariance of the potential V. Due to transla- 
tional invariance, the momentum ^ i is also conserved. 
Furthermore, the energy fluctuations are bounded due to 
the existence of a pseudo-Hamiltonian [l2T| . 



V. NUMERICAL TESTS 

A. Free rotation of an asymmetric body 

To demonstrate the advantage of the SIER integra- 
tors for an asymmetric body, simulations were per- 
formed on an isolated rigid body, whose moments of 
inertia are those of a water molecule, with initial con- 
ditions drawn from a canonical distribution at 297 K. 
For comparison, simulations were also performed us- 
ing the MRDL scheme of McLachlan, Reich, Dullweber 
and Leimkuhlerp], [§] and the so-called SEJ4 scheme of 
Celledoni and SafstromQ. The advantage of consider- 
ing free rigid body is that the exact motion is known 
and can be used to measure the deviation in the trajec- 
tories in addition to checking violations of energy con- 
servation. The energy is in fact insensitive to the ori- 
entation and the violations of energy conservation turn 
out to be small for all the integrations methods. To 
show that this does not mean that all methods are 
equally accurate, as an alternative error estimate, we 
used 6 = <±Tr[(A - A CX )(A - A cx ) T ]) 1 /2, with A ox the 
exact result 10]. While S m for A w A cx , 5 = 1 for a 
random rotation matrix A. Numerically it is found that 
the error scales as S = i/(10r) for moderate t, where r is 
the time at which the error in A becomes 10%. Hence r 
can be viewed as the time scale at which the orientation 
of the body is incorrectly calculated. This time scale 
depends on the time step h and the scheme used. For 
h = 1.66 fs, the time scale at which the orientation be- 
comes inaccurate is t « 117 ps using MRDL, while SEJ4 
gives r « 0.91 fis, and SIER2 gives r w 16 ms (SIER2 
and SIER4 have similar accuracy here). Given that in 
gases at standard conditions, typical free flight times are 
on the order of 100 ps, one might still consider the error 
in MRDL acceptable in such applications. However, time 
steps as large as 8 fs are possible for such systems [3] • For 
this time step, the time scale at which the orientation 
becomes inaccurate in the MRDL method is r pa 7.6 ps. 
Thus, the orientation would incorrectly given by MRDL 
well before a free flight is over, even though the violations 
in the energy conservation are small (~ 0.01%). Some- 
what surprisingly, then, the time step that can be used 
in MRDL simulations of a low density system is limited 
not by interaction parameters but by the accuracy of the 
free flight. On the other hand, SEJ4 gives r » 2.1 ns, and 
SIER2 gives r « 48 ms, so both give orientations which 
are still accurate after a free flight. In terms of computa- 
tional load, for a fixed time step, SIER2 was between 70% 
and 100% slower than MRDL, while SEJ4 was another 
20% slower still. Thus SEJ4 is slower and less accurate 
than SIER2. And while MRDL is roughly twice faster 
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h/f 




AH/AV in % 






AL Z /(L Z ) 




(fs) 


MRDL 


SEJ4 


SIER2 


SIER4 


MRDL(10 -13 ) 


SEJ4(10" 3 ) 


SIER2(KT 13 ) 


SIER4(10~ 13 ) 


4.2 


1.8 ±0.3 


0.088 ± 0.006 


0.088 ± 0.006 


0.0035 ± 0.0005 


16 ± 1 


0.0014 ± 0.0003 


6± 1 


2.6 ±0.7 


11 


11.4 ± 1.6 


0.55 ± 0.04 


0.55 ±0.4 


0.27 ±0.06 


7.0 ± 0.6 


0.069 ±0.019 


2.1 ± 0.7 


1.8 ±0.3 


21 


51 ±7 


2.24 ±0.16 


2.22 ±0.15 


14 ±2 


3.9 ±0.4 


2.3 ± 1.1 


1.4 ±0.3 


2.6 ± 0.4 



TABLE I: Conservation of 7i and L z for a dipole (1.84 Debye) in an electric field (2.7 MV/m), using various integrators. 



for given h, this cannot make up for the loss in accuracy. 



B. Water molecule in an electric field 

To test the SIER integrators on a non-free body, a 
constant electric field of 2.7 MV/m directed along the z- 
axis is introduced with which the water molecule interact 
via its dipole moment of 1.84 Debye. The strength of the 
electric field corresponds to that induced by a second wa- 
ter molecule at a distance of 7 A. Since the exact motion 
of the molecule is not known for this case, more conven- 
tional measures for accuracy are used. In TableJH the rms 
fluctuations of two conserved quantities, the total energy 
TL and the ^-component of the angular momentum vector, 
L z , are given for different time steps for the SIER inte- 
grators and compared with the MRDL scheme and the 
SEJ4 scheme. Note that the fluctuations of 7i are given 
relative to the fluctuations of the potential energy AV 
(a natural scale of energy fluctuations), while the fluc- 
tuations of L z are given relative to the average L z . For 
an equitable comparison, the relative accuracies of the 
simulations were compared as a function of h/f, where 
/ is the number of force evaluations per time-step. For 
MRDL, SEJ4, and SIER2, / = 1, while / = 4 for SIER4. 
We note that the relative real-time computational loads 
are comparable to those in the free case. 

It is evident from Table Q] that SIER2 has the same 
degree of energy conservation as SEJ4, outperforming 
MRDL in this respect. However, SIER2 and MRDL 
conserve L z equally well since the MRDL scheme also 
conserves the angular momentum exactly. On the other 
hand, the conservation of L z in the SEJ4 scheme is orders 
of magnitude worse than the other propagation schemes, 
casting doubt on the accuracy of its rotational dynamics. 
Thus, SIER2 combines the excellent energy conservation 
of SEJ4 with the exact angular momentum conservation 
of MRDL. In addition, SIER4 improves the energy con- 
servation for h/f below a threshold value of about 12 fs. 



C. Liquid water 

In simulations of high density liquids, one expects the 
free flight to be less important, and the advantages of us- 
ing SIER less. However, one still expects SIER4 to give 
better results than SIER2. To test this, simulations were 
also performed on a system of 512 water molecules at 
liquid density at a temperature of 297 K[20|. To ensure 



strict energy conservation, a smooth molecular cut-off 
of 11 A was employed. For the simulations, time steps 
ranging from h = 0.1 fs to h = 7fs were used (above 
which one sees a substantial energy drift). The results 
for SIER2 and SIER4 are shown in Fig.[TJ MRDL and 
SEJ4 simulations were also performed, but their results 
as well as their real-time computational loads were vir- 
tually indistinguishable from those of SIER2, and there- 
fore not plotted in Fig.[T] Thus, from the perspective of 
energy conservation, SIER2 performs as well as MRDL 
and SEJ4 here, although the trajectories are likely to be 
more accurate (as in the free case). From Fig.[TJ it is 
evident that below h/f « 1.8 fs, SIER4 is more accurate 
than SIER2, although its theoretical scaling of 0(h 4 ) for 
small h is frustrated by round-off errors 21 1 . 



VI. CONCLUSIONS 

In this paper, two symplectic, time-reversible algo- 
rithms for simulating rigid body dynamics were pre- 
sented, SIER2 and SIER4, which are of second and fourth 
order, respectively. The schemes do not use constraint 
conditions, nor an extended state space, nor Euler angles 
or quaternions, but instead make use of a recent imple- 
mentation of the exact solution of free rigid body mo- 
tion. The integrators conserve the symplectic structure 
on the conventional phase space of rigid systems and re- 
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FIG. 1: The rms energy fluctuations using SIER2 and SIER4 
for liquid water. All data are averaged over nine 1-ps runs 
divided by the fluctuations in the potential energy AV. 
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exactly when applicable. 

Numerical comparisons with two existing integration 
methods (MRDL and SEJ4) were performed for three 
systems: a free asymmetric rigid body, a water molecule 
in an external field, and liquid water. In the free case, a 
comparison with the exact trajectory was possible, which 
showed that the orientational dynamics given by the 
SIER schemes is superior to existing integration schemes 
at a given numerical cost. For the water molecule in 
an external field, the accuracy of the trajectories was 
assessed using the conserved quantities of energy and 
the z-component of the angular momentum. Here too, 
the SIER schemes performed better at a given numeri- 
cal cost. For the simulation of water at liquid densities, 
the energy conservation of the SIER methods is equal 
to that of the other integrators, due to the dominance 



of the forces in the accumulated error. Nonetheless, the 
trajectories are likely to be more accurate, because the 
free motion step in the integrator is performed exactly. 

The integrators presented here should be of use in as- 
trophysical simulations in which accuracy is important as 
well as in simulations of biomolecular and other systems 
in chemical physics, which often require propagation of 
many degrees of freedom to long time scales. 
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